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1  Problems  Studied 

The  primary  focus  of  this  project  was  on  the  analysis  and  development  of  new  parallel  algorithms 
for  the  solution  of  linear  and  nonlinear  initial/boundary  value  problems  (IB VPs)  in  two  space 
variables.  Parabolic,  second  order  hyperbolic,  biharmonic,  and  Schrodinger-type  problems  were 
considered.  The  new  algorithms,  which  are  alternating  direction  implicit  (ADI)  orthogonal  spline 
collocation  (OSC)  methods  employing  C1  piecewise  polynomial  spaces  of  arbitrary  order,  have 
been  implemented  and  their  efficacy  was  demonstrated  on  test  problems  taken  from  the  literature. 
Rigorous  stability  and  convergence  analyses  of  the  methods  were  also  carried  out. 

2  Summary  of  Results 

2.1  Parabolic  IB  VPs 

In  [3],  an  ADI  Crank-Nicolson  scheme  is  considered  for  the  solution  of  the  linear  parabolic  initial¬ 
boundary  value  problem 

du 

-x-  +  {Li  +  L2)u  =  f{x,y,t),  ( x,y,t)eClT , 

ot 

u(x,y,0)  =  gi(x,y),  (x,y)eCl, 
u(x,y,t)  =  g2{x,y,t),  (x,y,t)  e  dClr, 

where,  here  and  in  what  follows, 

Cl  =  (0, 1)  x  (0, 1),  ClT  =  x  (0, T],  dClT  =  dCl  x  (0, T], 
and  the  linear  differential  operators  L\  and  L2  are  given  by 

(1)  L\u  =  -a\ (x, +  h (x, y,  +  c(x,  y,  t)u,  L2u  =  -a2(x, y,t)^  +  b2(x, y,  t)^. 

In  comparison  to  the  scheme  outlined  in  the  original  proposal,  the  new  ADI  scheme  uses  in 

place  of  L2  and  L1^1  ■  Using  a  new  approach,  we  show  that  the  scheme  is  second-order  accurate  in 
time  and  of  optimal  third-order  accuracy  in  space  in  the  Hl  norm.  For  simplicity,  the  analysis  in  [3] 
is  presented  for  the  case  of  a  spatial  discretization  based  on  piecewise  Hermite  bicubics,  but  is  easily 
extended  to  OSC  discretization  with  piecewise  polynomials  of  higher  degree.  We  also  give  a  new 
efficient  implementation  of  the  scheme  and  test  it  on  a  sample  problem  for  accuracy  and  convergence 
rates  in  various  norms.  Earlier  implementations  of  ADI  OSC  schemes  were  based  on  determining, 
at  each  time  level,  a  two-dimensional  approximation  defined  on  Cl.  In  the  new  implementation,  at 
each  time  level,  we  determine  one-dimensional  approximations  along  horizontal  and  vertical  lines 
passing  through  Gauss  points  and  obtain  the  two-dimensional  approximation  on  Cl  at  the  final 
time  level  corresponding  to  t  =  T.  It  should  be  noted  that  with  respect  to  the  implementation  of 
OSC  schemes  the  non-divergence  forms  of  L\  and  L2  in  (1)  are  more  natural  than  the  divergence 
forms  of  L\  and  L2  which  are  typically  used  in  finite  element  (FE)  spatial  discretization.  In  fact, 
ADI  FE  Galerkin  methods  for  solving  variable  coefficient  parabolic  problems  in  the  divergence 
form  were  considered  in  [7],  However,  our  two-level,  parameter  free  ADI  OSC  scheme  does  not 
have  a  FE  Galerkin  counterpart.  The  method  of  [7]  of  comparable  accuracy  is  the  three  level  ADI 
Laplace-modified  scheme  requiring  the  selection  of  a  stability  parameter.  Our  ADI  OSC  scheme 
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with  piecewise  polynomials  of  degree  >  3  is  more  accurate  than  the  standard  ADI  finite  difference 
scheme  which  is  only  second-order  accurate  in  both  time  and  space. 

In  [4],  we  consider  a  nonlinear  parabolic  initial-boundary  value  problem  on  a  rectangular  poly¬ 
gon  with  the  solution  satisfying  variable  coefficient  Robin’s  boundary  conditions.  An  approximation 
to  the  solution  at  a  desired  time  value  is  obtained  using  an  alternating-direction  implicit  extrapo¬ 
lated  Crank-Nicolson  scheme  in  which  orthogonal  spline  collocation  with  piecewise  polynomials  of 
an  arbitrary  degree  >  3  is  used  for  spatial  discretization.  For  rectangular  and  L  shaped  regions  we 
describe  an  efficient  S-spline  implementation  of  the  scheme  and  present  numerical  results  demon¬ 
strating  the  accuracy  and  convergence  rates  in  various  norms.  For  problems  with  homogeneous 
Dirichlet  boundary  conditions,  we  observe  a  superconvergence  phenomenon  when  the  initial  con¬ 
dition  is  approximated  using  the  Gauss  interpolant  rather  than  the  quasi-interpolant  suggested  in 
[8]  for  parabolic  equations  in  a  single  space  variable. 

In  the  special  case  in  which  the  region  is  a  square  and  Dirichlet  boundary  conditions  are  pre¬ 
scribed,  the  problem  considered  in  [4]  is  of  the  form 


ai(x,y,t,u,  V«)- 


L7  U 

a2{x,y,t,u,Vu)-^>  =  f(x,y,t,u,Vu),  ( x,y,t )  G  fir, 


u{x,y,  0)  =  gi  {x,  y) ,  (x,  y)  G  ft, 


u{x,y,t)  =  g2{x,y,t),  (x,y,t)  G  dClT- 

We  have  carried  out  convergence  analysis  for  the  case  in  which  the  right  hand  side  /  depends  on  u 
and  Vu  but  the  coefficients  oq  and  a2  are  independent  of  these  quantities.  We  expect  to  extend  our 
analysis  to  the  case  in  which  a\  and  a2  depend  on  u  using  new  results  obtained  by  our  graduate 
student  Abdulrakhim  Aitbayev  who  in  his  Ph.D.  dissertation  [1]  obtained  new  convergence  results 
for  the  OSC  solution  of  the  nonlinear  elliptic  boundary  value  problem 


ai(x,y,u,  +  ai2(x,y,u,Vu)-^^  +  a2{x,y,u,Vu)^  +  c(a :,y,u,Vu)  =  f(x,y), 


(x,  y)  G 


u(x,y)  =  0,  {x,y)€d Cl. 


2.2  Hyperbolic  IBVPs 

In  [14],  two  schemes  are  formulated  and  analyzed  for  the  approximate  solution  of  the  linear  second 
order  hyperbolic  problem 

—-Z-  +Lu  =  f(x,y,t),  {x,y,t)e^r, 


u{x,y,0)  =  g0(x,y),  —(x,y,0)=gi{x,y),  (z,y)Gft, 

u{x,y,t)  =g2(x,y,t),  (x,y,t)  G  dClT, 


where  the  linear  differential  operator  L  is  given  by 


_  d  (  .  ^  <9u\  d  (  .  .cHA  .  .  ±.du  .  ,  du  .  . 

Lu  =  ~ di  r1 1 (x’ ' y’ di )  ~  (,°2 ^ y^fy)+hl (x’ V't)fc+  h ^ y,t^+  c(x’ y ’  )u- 
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OSC  with  piecewise  Hermite  bicubics  is  used  for  the  spatial  discretization.  The  resulting  system  of 
ordinary  differential  equations  in  the  time  variable  is  discretized  using  perturbations  of  standard 
finite  difference  approximation  procedures  to  produce  Laplace-modified  (LM)  and  ADI  schemes. 
It  is  shown  that  the  OSC  LM  and  ADI  schemes  are  unconditionally  stable  and  of  optimal  order 
accuracy  in  the  Hl  and  discrete  maximum  norms  for  the  space  and  time  variables,  respectively, 
provided  a  stability  parameter,  which  appears  in  all  LM  type  methods,  is  chosen  appropriately. 
In  the  convergence  analysis,  a  bicubic  interpolant  of  the  solution  u  of  (2)  is  used  as  a  comparison 
function.  The  algebraic  problems  to  which  these  schemes  lead  are  also  described  and  numerical 
results  are  presented  for  an  implementation  of  the  OSC  ADI  scheme  to  demonstrate  the  accuracy 
and  rate  of  convergence  of  the  method. 

In  [5],  we  consider  the  approximate  solution  of  the  linear  second  order  hyperbolic  problem 

-QP  +  (Li  +  L2)u  =  f{x,y,t),  ( x,y,t)enT , 

3u 

u{x,y,0)  =  g0{x,y),  —  {x,y,0)  =  gi(x,y),  (x,y)e£l, 

u(x,  y ,  t)  =  g2(x ,  y,  t),  (x,  y,  t )  €  dO.T, 


where  the  linear  differential  operators  L\  and  L2  are  given  by  (1).  With  Q  the  set  of  Gauss  points 
in  fl,  r  =  T/M,  and  tn  =  nr,  the  new  OSC  ADI  scheme  consists  in  finding  piecewise  Hermite 
bicubics  u%,  n  =  2, . . . ,  M,  such  that  for  n  =  1, . . . ,  M  —  1, 


(3) 


(/  +  0.5t2£IK({)  =  /({,(„)  +  2t-2<(0,  { e  g, 

(I  +  0.5rJ£5)K+1  +  <-')({)  =  {  €  S, 


where  u°h,  u^,  u% |an,  n  =  2, . . . ,  M,  are  assumed  to  be  given,  and  where  L",  are  the  differential 
operators  of  (1)  with  t  =  tn.  In  the  first  equation  of  (3),  for  each  Gauss  point  in  the  y-direction, 
n  —  1, . . .  ,M  —  1,  is  a  piecewise  Hermite  cubic  in  the  re-direction  satisfying 

fi»(a,e)  =  r-2(/  +  0.5T2£5)(u;+1+<-1)(o,O,  «  =  0,1. 


We  show  in  [5]  that  the  scheme  (3)  is  second  order  accurate  in  time  and  of  third-order  accuracy  in 
space  in  the  Hl  norm.  An  efficient  implementation  of  the  scheme  is  similar  to  that  for  the  scheme 
of  [3]  and  it  involves  representing  in  terms  of  basis  functions  with  respect  to  y  alone  while  u'l 
is  represented  in  terms  of  basis  functions  with  respect  to  x  only.  It  is  interesting  to  note  that  for 
variable  coefficient  hyperbolic  problems  our  parameter  free  ADI  OSC  scheme  (3)  does  not  have  a 
finite  element  Galerkin  counterpart. 


2.3  Biharmonic  Problems 

Our  graduate  student,  Lou  Zhuoming,  completed  his  Ph.D.  dissertation  and  graduated  in  May  1996. 
His  research  was  concerned  with  the  derivation  of  existence,  uniqueness  and  convergence  results 
for  OSC  methods  and  the  implementation  of  these  methods  for  the  solution  of  three  biharmonic 
problems,  based  on  the  splitting  principle  [15].  The  first  problem  comprises  the  biharmonic  equation 
in  f 1  with  u  =  g\  and  A u  =  g2  on  dfl.  This  problem  becomes  one  of  solving  two  nonhomogeneous 
Dirichlet  problems  for  Poisson’s  equation.  The  resulting  linear  systems  can  be  solved  effectively 
with  cost  0(N 2  log2  N )  using  the  matrix  decomposition  algorithm  of  [6].  In  this  case,  optimal  order 
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ift-norm  error  estates,  *  -  0,1,2,  are  derived  In 

on  the  horizontal  sides  of  80,  Au  to “  *eP  “  osc  Fourier  method  is  formulated  for  the 

^  cos,  0(W  loS2JV).  The  third  problem  -  the 

biharmonic  Dirichlet  problem 

a2«  =  f{x,y),  (x,y)en, 
u  =  gi{x,y),  ^  =  92(x,y),  ( x,y)eda , 

and  again  optimal  Hl-  and  H2- norm  error  estimates  we  denved  The  01 SC 

by  a  direct  method  which  is  based  <*J**^™\  ^s  capachanTe  matrix  algorithm  is 
problem  as  the  auxiliary  problem.  The  to  a  whkhP  in  partiCular,  verify  the  fourth  order 

Results  of  some  numerical  of  the  derivative  approximations  at  the 

accuracy  of  the  approximations  a  P  different  tools  to  simplify  the  error 

mesh  points.  The  paper  (16)  »  based  on  th rwort “ d  ^  questions. 

analyses  and,  in  addition,  answers  “me  outstandmg^  ^  ^  ^  December  1997.  His 

David  Knudson  completed  his  Ph.D.  d  *  1  Ciarlet-Raviart  mixed  finite  element 

research  was  concerned  with  the  piecewise  Hermffe cubic  C  „eous  boundary 

Galerkin  method  for  the  solution  of  the  bthannomc  D'-cWet  Ptoble  ^  osc  apptoach 

conditions.  The  aim  of  this  work  was  to  corner  ^ ‘he  °alerkm  ap Pf  ^  osc  methods 

considered  by  Lou  [15]  and  provi  e  insig  1  ,  ,  ^  existence  and  uniqueness  of  the 

for  this  problem  and  also  certain  Schrodinger  systems.  In  [9]  first  exist  Galerkin 

Galerkin  solution  were  proved.  Then  a  Schur  two  vertical  sides 

problem  to  a  Schur  complement  system  mvo  vi  g  PP  problem  with  An  instead  of 

of  and  to  an  auxiliary  Galerkin  problem  for  system  with  a  symmetric 

.  du/dn  specified  on  the  two  vertical  sides  of  «l.  The  Schur  caav^  J  t  method.  A 

and  positive  definite  matrix  was  solJe  j181^  ^oblem  for  a  related  biharmonic  problem  with  An 
preconditioner  was  obtained  from  the  Gale  P  conjectured  that  the  precondi- 

instead  of  Su/dn  specified  on  the  two  horizontal  s  des  o  »  “  “  ^  *  x  N  partili„n  the  cost 

tioner  was  spectrally  equivalent  to  the  Sc  ur  comp  e  ,he  schur  complement  matrix  by 

of  solving  the  preconditioned  system  and  the  cost  of  b  of  solving 

a  vector  are  0(«*)  each.  " ^  ^'^Won  to  the  auxiliSy  ialerkin  problem  is 
the  Schur  complemen  system  0(  J2^^  transforms  at  a  cost  of  0(N2  log2  N)  Hence 

obtained  using  separation  of  vana  w  n( jy2  logo  N).  Numerical  results 

~ st  iTandlu  L  Of  optimal  fourth 

indicate  that  the  L  and  H  norm  er  vv  f  ,  d  for  t^e  approximations  to  u 

and  third  orders  respectively.  Convergence  a,  the mote ^  fourth  ord  ^  PP  Au  A 

and  An  and  third  order  for  the  approximations  to  the  first  order  denvau 

based  on  this  research  is  being  prepared  for  publication  .  bicubic  OSC  solution  of  the 

In  [16],  the  Unear  system  °a  ca  acltance  matrix  algorithm  whose  cost  is  0(N3). 

biharmonic  Dirichlet  problem  was  so  g  developed  for  solving  this  linear  system. 

linear  system  using  the  preconditioned  conjugate  gradient  metho  . 
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2.4  Schrodinger  Systems 

In  [11],  Crank-Nicolson  and  ADI  OSC  schemes  are  formulated  and  analyzed  for  the  approximate 
solution  of  the  linear  Schrodinger  problem 

^7  -  iAip  +  ia(x,y,t)ip  =  f(x,y,t ),  (x,y,t)  € 

at 

(4)  ip(x,  y,  t)  =  0,  (x,  y,  t)  €  dtir, 

ip(x,y,0)  =  ip°{x,y),  (x,y)€tl, 

where  i2  =  -1,  A  is  the  Laplacian  and  a  is  a  prescribed,  real  function,  while  ip,  ip°  and  /  are 
complex- valued.  This  type  of  problem  arises  in  many  disciplines,  such  as  quantum  mechanics, 
underwater  acoustics,  plasma  physics  and  seismology.  Usually  /  =  0  in  quantum  mechanics.  We 
write  ip,  f  and  ip°  as  ipi  +  iip2,  fi  +  if2  and  ip\  +  iip%,  respectively.  Taking  real  and  imaginary  parts 
of  (4)  then  yields 

du 

—  +  S(-A  +  a(x,  y,  t))u  =  F(x,  y,  t),  (x,  y,  t )  €  fir, 

(5)  u(x,y,t)  =  0,  ( x,y,t )  €  8Qt, 


u(x,y,0)  =  u°(x,y),  (x,y)  G  0, 


where 


u 


=  (  i>i 
1  ^2 


,  F  = 


and  V?=(  ^ 


/2 


r2 


are  real-valued  vector  functions,  and  S  is  the  2  x  2  skew-symmetric  matrix 


0  - 
1  0 


Hence,  (5),  and  thus  (4),  is  not  parabolic  but  a  Schrodinger-type  system  of  partial  differential 
equations.  In  this  paper,  OSC  with  Cl  piecewise  polynomials  of  arbitrary  degree  r  >  3  in  each  space 
variable  is  used  for  the  spatial  discretization  of  (5).  The  resulting  system  of  ordinary  differential 
equations  in  the  time  variable  is  discretized  using  the  trapezoidal  rule  to  produce  the  Crank- 
Nicolson  OSC  scheme,  which  is  then  perturbed  to  obtain  the  ADI  OSC  scheme.  The  stability  of 
these  schemes  is  examined  and  optimal  order  a  priori  error  estimates  in  both  the  H 1  -norm  and 
the  L2-norm  at  each  time  step  derived.  The  approach  employed  in  the  convergence  analysis  of 
the  schemes  is  based  on  using,  as  a  comparison  function,  a  projection  of  the  exact  solution  into  a 
space  of  C 1  piecewise  polynomials  of  degree  r  in  each  space  variable.  This  is  a  key  element  in  the 
derivation  of  the  optimal  order  L2-error  estimates. 

In  [12,  13],  OSC  methods  are  considered  for  the  solution  of  problems  governed  by  the  equation 


(6) 


+  2^  +  A  {q{x,  y)Au)  =  f(x,  y,  t), 


(x,y,t)  €  Ot, 


where  v  is  a  nonnegative  constant  related  to  viscous  damping  coefficient  (u  =  0  if  external  viscous 
damping  is  ignored),  and  q(x,  y)  is  a  variable  density  function  such  that 


0  <  9min  ^  V)  —  9max  <  00>  {xil/)  £ 
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see  [10].  This  equation  commonly  arises  in  plate  vibration  and  seismological  problems.  The  initial 
conditions  are 

u(x,y,0)  =  go{x,y),  0)  =  gi(x,y),  (x,y)eSl, 

while  the  boundary  conditions  (BCs)  are  of  one  of  the  following  three  types: 

BC1:  “clamped”  BCs: 

du 

u(x,y,t)  =  0,  —  (x,y,t)-0,  (x,  y,  t)  G  d9.T\ 

BC2:  “hinged”  BCs: 

u(x,  y,  t)  =  0,  Au(x,  y,t)  =  0,  (x,y,t)  €  dtlr', 

BC3:  BCs  in  which  the  vertical  sides  are  hinged  and  the  horizontal  sides  are  clamped: 

u(x,y,t)  =  0,  (x,y,t)  G  dClT, 

A u(x,y,t)  =  0,  {x,y,t)  G  5Qi  x  (0,T], 

du 

—  (x,y,t)  =  0,  (x,y,t)  GdD2  x  (0,T], 

where 

dft\  =  {(a,  y)  :  a  =  0, 1,  0  <  y  <  1},  dti?  =  {(x,  a)  :  0  <  x  <  1,  a  =  0, 1}. 

We  reformulate  each  problem  by  introducing  the  functions 


For  example,  setting 


Ul 


U2  =  qeutAu. 


U  = 


G  = 


9i 

qAg0  ’ 


in  the  problem  with  BC1,  we  obtain  the  Schrodinger  system 


where 


^  -  SgAU  +  vRU  =  F,  (x,  y,  t)  G 

<  U(x,y,0)  =  G{x,y),  (x,y)efl, 

3v»i 

ui(x,y,t)  =  0,  ~(x,y,t)  =  0,  {x,y,t)  G  <9fir, 


S<  = 


0 

q 


-l 

o  ’ 


R  = 


1  0 
0  -1 


Note  that  it  is  necessary  to  carry  out  an  additional  calculation  to  obtain  an  approximation  to  u 
from  that  to  u\. 

To  determine  an  approximation  to  U,  we  used  OSC  with  C 1  piecewise  polynomials  of  arbitrary 
degree  r  >  3  in  each  space  variable  for  the  spatial  discretization.  The  resulting  systems  of  ordinary 
differential  equations  in  the  time  variable  were  then  discretized  using  standard  Crank-Nicolson  or 
ADI  techniques  involving  only  two  time  levels  (cf.  [11]).  Specifically,  we  formulated  and  analyzed 
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Crank-N  .colson  OSC  schemes  for  ah  three  choices  of  BCs  mrd  £  “^“"1 

for  BC2  and  BC3  and  a^lyaed  “"it  Dixichlet  problem  (16,  we 

v  =  0.  Based  on  our  experience  with  OSC  metn  method  for  BC1.  For  this  problem, 

believe  that  it  is  not  possible  to  formu  ate  a  s  an  Crank-Nicolson  scheme 

we  employed  the  capacitance  matrix  method  (d.  [If , polynomials  on  a  unif0rm 
efficiently  for  the  the  special  case  o  cons  an  d  stability  of  each  scheme  and  show  that, 

partition  of  a  We  is  of  optimal  order  r  +  1  -  m  in 

“d‘l“  n^e.’  L^ntationa,  issues  were  afco  addressed  and  nnmerrcal  results 

obtained  which  confirm  the  theoretical  analyses. 
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